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ABSTRACT 

Context. The migration of planets plays an important role in the early planet-formation process. An important problem has been that 
standard migration theories predict very rapid inward migration, which poses problems for population synthesis models. However, it 
has been shown recently that low-mass planets (20-30 M^o,,;,) that are still embedded in the protoplanetary disc can migrate outwards 
under certain conditions. Simulations have been performed mostly for planets at given radii for a particular disc model. 
Aims. Here, we plan to extend previous work and consider different masses of the disc to quantify the influence of the physical disc 
conditions on planetary migration. The migration behaviour of the planets will be analysed for a variety of positions in the disc. 
Methods. We perform three-dimensional (3D) radiation hydrodynamical simulations of embedded planets in protoplanetary discs. 
We use the explicit-implicit 3D hydrodynamical code NIRVANA that includes full tensor viscosity, and implicit radiation transport. For 
planets on circular orbits at various locations we measure the radial dependence of the torques for three different planetary masses. 
Results. For all considered planet masses (20-30 Me„„i,) in this study we find outward migration within a limited radial range of 
the disc, typically from about 0.5 up to 1.5-2.5 aj„,,. Inside and outside this interval, migration is inward and given by the Lindblad 
value for large radii. Interestingly, the fastest outward migration occurs at a radius of about aj„,, for different disc and planet masses. 
Because outward migration stops at a certain location in the disc, there exists a zero-torque distance for planetary embryos, where 
they can continue to grow without moving too fast. For higher disc masses (Mf;,„, > 0.02Mo) convection ensues, which changes the 
structure of the disc and therefore the torque on the planet as well. 

Conclusions. Outward migration stops at different points in the disc for different planetary masses, resulting in a quite extended 
region where the formation of larger cores might be easier In higher mass discs, convection changes the disc's structure resulting in 
fluctuations in the surface density, which influence the torque acting on the planet, and therefore its migration rate. Because convection 
is a 3D effect, 2D simulations of massive discs with embedded planets should be handled with care. 

Key words, accretion discs - planet formation - hydrodynamics - radiative transport - planet disc interactions 



1. Introduction 

Planetary migration of embedded low-mass planets in isother- 
mal discs indicates inward migration, so that the planet 
might be los t in th e star before the accretion disc is gone 
dTanaka et alj 2002[). Recent studie s (starting with the work 
of iPaardekooper & Mellemal ( l2006l) ) have shown that the in- 
clusion of radiation transport in planet-disc interaction stud- 
ies resulted in a slowed down or even reversed migration 
dBaru teau & Massetl l2 008al: IPaardekooper & Papaloizoul 120081: 
Paardekooper & Melle mall2008t iKlev & Cridall2008t iKlev et all 
200^iAvliffe & Batai2oTon20m 

All authors agree that the inclusion of radiation transport is 
an important effect that should be considered, however, not all 
authors find outward migration. The efficiency of outward mi- 
gration depends on the magnitude of positive corotation torques. 
These are determined by the local entropy and vortensity gradi- 
ents. For isothermal and adiabatic discs these gradients cannot 
be maintained and so-called torque saturation reduces the coro- 
tation effects. However, the inclusion of radiative transport and 
viscosity prevents saturation, and the negative Lindblad torques 
(caused by the spiral arms) can be overwhelmed by the posi- 
tive contributions from the corotation regio n. The magnitude of 
the effect depends on the planetary mass (iKlev & Cridall20()8t 
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iKlev et all 20091) and the tempe rature gradient in the disc. For 
example, lAvlifFe & Batd (1201 ll) found for temperature slopes 
of j6 > 1 .0 (with T oc r^^) that outward migration is possi- 
ble even for p lanets with up to SOMEanh in 3D simulations. 
The migration rate increases with an increasing temperature 
slope. T hese results are also refl ected by the theoretical analy- 
sis from 'Ma sset & Casolil (|201(A. Recently, improved theoreti- 
cal torque formulae for low-mass planets embedded in an adia- 
batic disc have been prese nted bv iMasset & Casolil (l2009h and 
IPaardekoope r et aP (l2010l) . The most recent improvements con- 
sider the nec essary inclusions of radiative diffusion and viscosity 
(fMasset & Casolil 120101: IPaardekooper et al.ll20Tlh . These have 
been developed for 2D discs where the diffusive effects can oper- 
ate only in the disc's plane. Interestingly, radiative cooling alo ne 
can lead to similar, even stronger effects dKlev & Cridall2008h . 

A very important open question is how far the planet will 
move outwards in a fully radiative disc. When the protoplanets 
are stopped from migrating outwards at a certain point in the 
disc, a protoplanet moving inwards from farther out will stop 
at the same location in the disc. Therefore, a planetary trap in- 
side the disc can be created, which can act as an area of plane- 
tary mergers, leading to larger cores. A planetary trap inside the 
disc created by surface density changes can function as a feeding 
zone to planetary cores (Morbidelli et al. 2008), but it is unclear 
how realistic these surface density changes are. The inclusion of 
radiation transport/cooling in a disc might provide such a trap in 
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a normal disc structure, at least for planets within a given mass 
range. ^ 

In ISandor et alj ( 1201 lb the authors show in N-body sim- 
ulations th at because of the inclus ion of unsaturated type-I- 
migration jPaardekooper et al]|2010l) large planetary cores (up 
to lOMEarth) Can form in protoplanetary discs well before the 
disc is accreted onto the star 

In population synthesis models that include the standard mi- 
gration prescription for isothermal discs the migration is too fast, 
such that the outcome distribution does not match the observed 
one. Specifically, the type-I-migration rate should be 10 to 1000 
times slower than exp e cted from linear analy sis jAlibert et al.l 
l2004t llda&Linll2008l: lMordasiniet"ani2009l) . As pointed out 
above, this could be provided by the inclusion of radiation 
transport in hydrodynamical type-I-migration simulations. An 
upda ted version of the type-I-m igration, by using the formula 
from lPaardekooper et al.l (1201 Cj). in population synthe sis models 
shows very promising results dMordasini et al.ll2010h . but addi- 
tional simulations are needed to verify these a ssumptions. 

In our previous simulations and studies ( iKlev et alj 120091 : 
iBitsch & Kiev 2010; Bitsch & Kley 2011]) we have only con- 
sidered one standard disc model with a given mass. In reality, 
of course, protoplanetary discs can have a variety of masses. A 
fully radiative disc will settle in an equilibrium state between 
viscous heating and radiative transport/cooling. Given that all 
other physical parameters are fixed, the resulting disc structure 
only depends on the disc's mass and viscosity, because a more 
massive or more viscous disc is able to heat the disc more effec- 
tively. A more massive disc therefore influences the migration 
rate of embedded planets. Here we focus on low-mass planets 
that can undergo outward migration in fully radiative discs. 

In this paper we extend previous studies and investigate the 
possible radial range over which outward migration can occur 
and analyse the influence of the disc's mass on the migration in 
detail. An important effect is the onset of convection in the disc, 
which becomes stronger for more massive discs. 

In Section |2] we give a short overview of our co de and nu- 
meric al methods, where more details can be found in (iKlev et al.l 
|2009|) . The influences of the distance of the embedded planet to 
the central star is discussed in Section |3] We then analyse the 
influence of the disc's mass on the disc's structure (density, tem- 
perature, aspect ratio) and then the influence on migration of em- 
bedded low-mass planets in Section |4] Convection in the disc is 
also briefly discussed in Section|4] We then summarize and con- 
clude in Section|5] 

2. Setup of the simulations 

The protoplanetary disc is modelled as a three-dimensional 
(3D), non-self-gravitating gas whose motion is described by the 
Navier-Stokes equations. We treat the disc as a viscous medium, 
where the dissipative effects can then be described via the vis- 
cous stress-tensor approach. We also assume that the heating of 
the disc occurs solely through internal viscous dissipation and 
ignore the influence of additional energy sources (e.g. irradiation 
form the central star). This internally produced energy is then ra- 
diatively diffused through the disc and eventually emitted from 
its surface. For this process we use the flux-limited diff'usion ap- 
proximation, which allows us to treat the transition from opti- 
cally thick to thin regions as an approximation. A more detailed 
description of the modelling an d the numerical methodology is 
provided in our previous pa pers (iKlev e t al. 2009: B itsch & Kiev 
I2OIOI: iBitsch & Klevll201ll) . and for that purpose we limit our- 
selves here to present only the most necessary information. 



2.1. Physical setup 

We solve the Navier-Stokes equations numerically using a spa- 
tially second order finite volume method based on the code 
NIRVANA (Ziegler & Yorke 1997), with impUcit radiative trans- 
port in the flu x-limited diffusion approx imation and the F ARGO 
(lMasselll2000l) extension as described in lKlevetall (l2009l) . We 
use a spherical polar coordinate system (r, 6, <p), where the com- 
putational domain consists of a complete annulus (0 < < 2n) 
of the protoplanetary disc centred on the star, extending from 
r,„in to r„„v, which are determined by the location of the planet. 
For symmetry reasons and because we only use non-inclined 
planets, we restrict ourselves in the standard setup to the up- 
per half of the disc. Hence, in the vertical direction the annulus 
extends from the equator up to 7° above the disc's midplane, 
meaning 83° < < 90°. Here 9 denotes the polar angle of our 
spherical polar coordinate system measured from the polar axis. 
For the study of convection we use in addition a two-sided disc, 
see below. The central star has one solar mass M» = Mq, and 
the total disc mass inside [r,„,„,r,„av] is Mdisc - O.OIMq, un- 
less stated otherwise in Section |4] For our radiative models the 
temperature and subsequently Hjr is calculated self-consistently 
from the equilibrium structure given by the viscous internal heat- 
ing and radiative diffusion. We note that for given physics (equa- 
tion of state, opacity, viscosity) the structure of the disc is solely 
dependent on its mass, and this is one aspect that we will inves- 
tigate in this paper. 

For the radiative transport we use a one-temperature ap- 
proach and apply the flux-limited diffusion a pproximation usin g 
analytic Rosseland opacities, for details see iKlev et al] ( |2009|) . 
To close the basic system of equations we use an ideal gas equa- 
tion of state with constant mean molecular weight yU - 2.35 for 
a standard solar mixture. The adiabatic index is y = 1.4. For 
the present study, we use a constant kinematic viscosity coef- 
ficient with a value of v = 10'^ cm^/s, a value that relates to 
an equivalent a - 0.004 at 5.2AU for a disc aspect ratio of 
Hjr - 0.05, where v = aH^Q.K- In standard dimensionless units 
with ro - ajiip - 5.2AU and fo - ^^^'('"0) ' we have v = 10"^. 

2.2. Calculation setup 

Our coordinate system rotates at the initial orbital frequency of 
the planet (at r - rp). We use an equidistant spherical grid 
in r,0,4> with a resolution of (Nr,Ne,N^) = (266,32,768) ac- 
tive cells for our simulations. At r„„„ and r,„a_^ we use damping 
boundary conditions for all three velocity components to min- 
imize disturbances (wave reflections) at these boundaries. The 
velocities are relaxed towards their initial state on a time scale of 
approximately the local orbital period. The angular velocity is 
relaxed towards the Keplerian values, while the radial velocities 
at the inner and outer boundaries vanish. Reflecting boundary 
conditions are applied for the density and temperature in the ra- 
dial directions. We apply periodic boundary conditions for all 
variables in the azimuthal direction. In the vertical direction we 
set outflow boundary conditions for 0„„„ - 83° (the surface of 
the computational domain). 

In our previous work, we have discussed the calculation of 
the torque acting on a planet in great detail. Outward migration 
seems only possible and is strongest when the planet is on nearly 



seems only possible ana is strongest wnen tne planet is on nearly 
circular orbits in the midplane of the disc dBitsch & Kiev 20101: 
'Bitsch & Klevll201lh . Additionally, we stated in IBitsch & KlevI 
(2011) that the measured migration rate from planets on fixed 
orbits is equivalent to the migration rate determined from mov- 
ing planets in discs. Hence, we consider here primarily planets 



Bitsch & Kley: Range of outward migration and influence of the disc's mass on planetary migration 



on fixed circular orbits in the midplane of the disc, and calculate 
the torque acting on the planet, because the torque represents a 
direct measurement of migration in this case. For a one-sided 
disc we use symmetric boundaries at 6,„ax - 90° (the disc's mid- 
plane). To correctly see the influence of convection in the disc 
we use two-sided discs for some high-mass discs. For these sim- 
ulations we used outflow boundaries for both 0„„„ and 9„,ax- 

The models are initialized with constant temperatures on 
cylinders with a profile T{s) oc 5 ' with s - rsinft The ini- 
tial vertical density stratification is approximately given by a 
Gaussian where the radial surface density follows a Y.{r) oc r"'^^ 
profile. In the radial and ^-direction we set the initial velocities to 
zero, while for the azimuthal component the initial velocity u^ is 
given by the equilibrium of gravity, centrifugal acceleration, and 
the radial pressure gradient. This corresponds to the equilibrium 
configuration for a purely isothermal disc. Before embedding the 
planet, we bring the disc into radiative equilibrium by perform- 
ing first 2D axisymmetric simulations in the r — 9 plane. This 
takes about 100 orbits. We then extend this model to a full 3D 
simulation by expanding the grid into the (;i-direction, and the 
planet is embedded. 

For the gravitational potential of the planet we utilize the 
cubic potential, where the potential is exact beyond a transi- 
tion (smoot hing) radius rsm and smoothed by a cubic polyno- 
mial inside (iKlahr & Klevll2006HKlev et al.ll2009h . Here we use 
Vsm - O.SRh, where Rh is the Hill radius of the planet. 

The torques acting on 20, 25, and 30MEarih planets are calcu- 
lated to determine the direction of migration. The gravitational 
torques acting on the planet are calculated by integrating over the 
whole disc, where we apply a tapering function to exclude the in- 
ner parts of the Hill sphere of the planet (Crida et al. 2008) . This 
torque-cutoff is necessary to avoid strong, probably noisy contri- 
butions from the inner parts of the Roche lobe and to disregard 
material that is gravitationally bound to the planet dCrida et al.l 
|2009|) . Here we assume (as in our previous papers) a transition 
radius of 0.8 Hill radii. 

3. Range of outward migration 

Previo us simulati ons by several authors ( Baruteau & Massetj 



2008a: 



^ ■2008t 

Paardekooper & MeUemall2008l: iKlev & Cridall2008l: iKlev et al.l 



Paardekooper & Papaloizo 



20091) indicated that outward migration of low-mass planets 
may be possible during an early evolutionary state of planet 
formation. However, because the simulations dealt mostly with 
planets at a given distance from the star, typically 5.2 AU, the 
radial range over which the migration ma y be directed outwards 
has not been addressed in great detail. In lBitsch & Kiev! (1201 Ot) 
we obtained some results for moving planets but the extent of 
the outward migration remained unclear. 

In order to address this problem, we simulate 20,25, and 
3QMEarih plancts on fixed circular orbits embedded in fully ra- 
diative discs at various distances from the host star. The planet's 
semi-major axis rp lies in a range of 0.6 < rp < 5.0ro, where 
f() = ojup = 5.2AU. With increasing distance from the star, the 
density and temperature of the disc decrease and at some point 
in the disc the conditions for outward migration might not be 
fulfilled anymore. 

For this set of simulations with different planet locations we 
use a disc setup with a density profile such that the total disc 
mass equals Moisc = O.OIM© for a planet at rp - 1 and r,„/„ = 0.4 
and r„,„i - 2.5 in units of r^. The planets are embedded in a way 
that the distance to the inner edge is always the planets loca- 
tion divided by 2.5, while the distance to the outer edge is the 



planet's location to the star multiplied by 2.5. This setup ensures 
that the radial boundaries are always far enough away so they do 
not influence our results of embedded planets. Since the overall 
surface density profile (E oc r""^) of the different disc models 
refers to the same physical situation, the total disc mass in the 
computational domain changes in the same way as the computed 
domain changes its size. The surface density at a given distance 
to the central star is constant for all disc models. The rotation 
frequency of the grid matches the rotation speed of the planet, 
so that the planet remains at a fixed position inside the computa- 
tional grid at all times. 

In Fig. [1] the specific torque (per planet mass) acting on the 
three planetary masses at different distance from the central star 
is displayed. For all planet masses the most extended positive 
torque (indicating outward migration) is around r x 1.0fly„p. 
At longer distances to the central star, the torque acting on the 
planets decreases continuously to negative torques, and this tran- 
sition from positive to negative torques occurs at larger dis- 
tances for lower planet masses. For the lowest mass planet with 
IQMEarth the transition is at r a; l.Aajup (zero-torque distance 
to the central star). With even longer distances the torques re- 
main negative but with diminishing strength, indicating inward 
migration. For higher planetary masses (25 and SOM^a,,/,) the 
zero-torque distance is decreasing (1.9 and lAcijup). For shorter 
distances to the central star, inside the maximum, the torque act- 
ing on the planet is reduced again until it reaches about zero for 
rp - O.Sajup for all planetary masses. 

In Fig.|2]the torque for the 20M Earth plan et is again displayed 
together with semi-analytical results from [Paardekooper et aU 
(l20I0i l20IIh. The black boxes refer to the torque formula pre- 
sented in [Paardekooper et aP (l2010l) . which applies for the un- 
saturated torques in adiabatic discs. It is given by 

rr,o,/ro = -2.5 - 1.1/3 + O.la + 1.1(1.5 - a) + 7.9^/r , (1) 

with a and /? being the slope of the radial surface density and 
midplane temperature profile, respectively 



S(r) oc r ", T{r) oc r 



(2) 



To calculate the semi-analytic results in the plot we use a con- 
stant a = 0.5, and/? changes from about 1.7 at r = l.Ofly,,^ to 
0.33 at r = S.Oajup, see also Fig.Qbelow. The slope of the en- 
tropy profile is then given by ^ = /3-(j-l)a - 1 .5 at r = 1 .Qaj„p. 
The torque is normalized to 



(!/ 



ipr^nl 



with q the planet/star mass ratio, h the relative disc height (H/r), 
Ep the surface density at the planet's location and Q.p the rotation 
frequency of the planet in the disc. The vertical height H of the 
disk is defined as H - Cs/Q. where the midplane values are used 

for the adiabatic sound speed c, and the angular velocity il. 

At r = 1.0a jup the formula from [Paardekooper et al.l (1201 Oh 
agrees well with our 3D simulations. As can be inferred from 
Eq. ([1]), the theoretical torque could never become negative for 
constant slopes (a,/3,^). However, as /3 varies with r, a change 
from positive to negative torques is possible and indeed occurs 
at r X 3.5 (black squares in Fig. |2]i. However, in the range of 
1.2 < r < 3.8 the torque formula predicts a much higher torque 
than our 3D simulations, as is the case for larger distances to the 
central star, where the formula suggests a larger negative torque 
compared with our results. It should be noted that Eq. ([T]i in- 
cludes Lindblad and corotation torques where the latter include 
contributions from the vorticity as well as entropy gradient. One 
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Fig. 1. Torque acting on planets with three different masses em- 
bedded in fully radiative discs located at various distances from 
the star (red for IQMEariiu blue for 25 M Earth and purple for 
30MEanh)- The torques are plotted when planet and disc have 
reached a quasi-stationary equilibrium where the torque acting 
on the planet is approximately constant in time. 



Fig. 3. Torque acting on a lOMEanh planet in fully radiative discs 
as a function of distance from the star The additional curves in- 
dicate the esti mates from the theoretica l formula for viscous, dif- 
fu sive disks by Masset & Casolil (1201 Ol) and for isothermal disks 
bv lD'Angelo & Lubowl(l2010h . 
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Fig. 2. Torque acting on a IQMEanh planet in a fully radia- 
tive disc as a function of distance from the star (red plus 
signs). Additional c urves indicate the results from the torque 
for mula published in|Pa ardekoo peret al.l (1201 Ol) (black squares) 
and lPaardekooper et all (,201 1.) (purple cross es). Additionally we 
also p lot the Unear Lindblad torque by the iPaardekooper et all 
(1201 ih formula (blue asterisks). 



should also be aware that Eq. ([T]i is valid only at the beginning 
of the evolution when the torques acting on the planet are not 
saturated. However, at later times the flow settles to an equilib- 
rium state where the torques saturate and Eq. ([T]i is not valid any 
more. 

In Fi2.|2]we a dditionally plot re sults from the improved for- 
mula in IPaardekooper et al.l (1201 ll) (purple crosses). Its deriva- 
tion is quite complex and we give a brief summary in the 
appendix. The formula captures the behaviour of the torque 
caused by Lindblad resonances and horseshoe drag on low-mass 
planets embedded in gaseous discs in the prese nce of viscous 
and the rmal diffusion (iPaardekooperet al.ll201ll) . The new for- 
mula of IPaardekooper et al.l (1201 ll) includin g diffusion provides 
a slig htly better fit compared to the adiabatic lPaardekooperet aP 



a slig 
(I201C 



2010) formula. The torque is positive for r < 1.8ay„p and be- 
comes negative for larger distances to the central star In general 
the formula captures the trend of our simulations, but lies con- 



sistently below our radiative disks. Also the zero-torque equilib- 
rium radius differs. At the often used reference radius r - ajup 
the torque formula and our 3D simulations differ by a factor of 
about three. For r > 2.5aj„p our simulations show a negative 
torque acting on the planet, which continues to grow until about 
r s! A.Oajup. For these larger radii the difference between our 
results and the analytic formula increases. 

In Fig. |2l t he L indblad torque according to 
IPaardekooper et all ( 1201 lb is also displayed (blue). For 
larger distances to the star, the disc becomes thinner and the 
effects of the corotation torque become less important. The 
total torque should therefore converge towards the Lindblad 
torque, however the torques from our simulations differ by 
a factor of three at r - 5.Qaj,in with the Lindblad torque of 
IPaardekooper et all ( 1201 lb . As the Lindblad torque is dependent 
on the temperature gradient P, Eq. ([B, and on a pre- factor, a 
reduction of the pre-factor as in lMasset & Casolil (1201 Ol) reduces 
the Lindblad torque (see also Fig.jSj. 

In Fig. |3] we also plot the results from iMasset & Casolil 

(1201 oh . who provided an alternative recipe for planetary migra- 
tion in viscous discs with thermal diffusion. We do not cite the 
full formulae here, due to their complexity. In paragraph 8 of 
iMasset & Casolil (l2010l) . a summary of the torque formula is 
gi ven. Additionally, the isothermal results of the 3D simulations 
of lD'Angelo & Lubowl ( l2010h are plotted. 

Compared t o our 3D results the torque formula of 
IMasset & Casolil (l2010l) leads to qualitatively different results. 
The torque predicted by their formula is negative for all in- 
vestigate d distances to the central star. For large distances the 
torque of IMasset & Casolil (1201 Ol) matches better with our sim- 
ul ations, but is still about 35% too low. The Lindblad torque 
of IMasset & Casoli ( 20101) is smaller compared to the Lindblad 
torque in IPaardekooper et all (1201 Ol) due to a different pre- 
factor in front of the tem perature slope P (in Eg. ([p ). The 
3D isothermal formula by iD'Angelo & Lubowl (1201 Ol) shows 
a larger Lindbla d torqu e compared to our simulations and to 
IMasset & Casolil (l2010l) . But as the formula is for isothermal 
discs only, the torque has to be reduced by a factor of y compared 
to our radiative simulations or compared to the torque from 
IMasset & Casolil ( l2010l) . Taking this into account the torques 
match quite well for large distances to the central star. 
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To test our i mplem entation of the torque formula by 
Masset & Casoli (l2010l) . we applied it to the results of 
AvlifFe & Bate (20 111) using our temperature gradient ( see their 
Fig. 3). However, in contrast to lAvUffe & Bate! (1201 ih we find 
negative torques when applying that formula. Apparently, when 
comparing their numerical results to the anytical formula by 
iMasset & Casolil (l2010l) . Ayliffe & Bate applied the latter with a 
somewhat larger viscosity (approx. by a factor of In) than quoted 
in text. This resulted in the positive torque quoted in their paper. 
(B. AylifFe, private communication). 

There are many reasons for the differenc es between the sim- 
ulations and the for mulae. The formulae in Paarde kooper et aP 
(120101 1201 lb and in iMasset & Ca"soiil ( 12010) were derived for 
2D discs only, but the horseshoe drag can be eve n stronger in 
three d imensions, as shown for isothermal discs in IMasset et al.l 
(|2006|) . The radiative diffusion is also most effective in the verti- 
cal direction, meaning that the two-dimensional approximation 
where diffusion is is restricted to the o rbital plane, does no t cap- 
ture the true physical effects (see also lKlev & Cridall2008l) . The 
formulae were derived for background gradients of temperature 
and surface density, but as the disc with an embedded planet 
evolves, the structure of the disc changes and the basic assump- 
tions (gradients in temperature, density, and so on) used to derive 
the formulae might not be valid any more. This is in particular 
true for planets in the mass range studied here, because the the- 
ory has been developed for lower mass planets, which do not al- 
ter their surroundings significantly. A more detailed discussion 
about the smoothing of the planet in iPaardekooper et al] (1201 lb 
can be found in Appendix |A] Nevertheless, there seems to be a 
very rough qualitative agreement between the improved torque 
formula of Paardekooper et al. ( 2011) and our numerical results. 

In Fig. |4] the radial torque density r(r) is displayed for 
IQMEarih plancts at rp - 0.6, rp = 1.0, rp - 2.5 and rp - 
A.Oajiip, where the radius is normalized to the actual planetary 
distances to allow a direct comparison. In all curves the contri- 
butions by the Lindblad torques are clearly visible, positive for 
r < 1 and negative for r > 1 . The contribution responsible for 
the torque reversal, the 'spike' just inside r = 1, is visible only 
for the rp = 0.6 and rp - 1.0 locations, which both show an 
overall positive total torque. Even though the torque acting on 
the rp = Q.bajup planet is much reduced, the spike in the torque 
distribution is clearly visible. In lKlev et al.l ( l2009l) we discussed 
the origin of the spike for the rp - \ .Oajup planet. It is an indica- 
tor for a density enhancement ahead of the planet just inside of 
rp, and thus creating a positive torque. 

For the rp = 2.5a j^p case the corotation spike in the fully ra- 
diative case in the torque density is not visible any more, only the 
Lindblad torques are visible. The resulting total torque is about 
zero, which indicates that the Lindblad torques are just counter- 
balanced by the corotation effects. For even longer distances to 
the central star, the torque is identical to the (negative) Lindblad 
torque , indicating inward migration (see also Fig.[T]i- 

In iKlev et all (l2009l) we explained in great detail with the 
help of 2D surface density plots how the torque acting on 
the planet is created in fully radiative discs. The torque act- 
ing on a lOMEarth planet on a circular orbit at rp = ay,,^ in 
a fully radiative disc is positive, indicating outward migration. 
In Fig. |5] the 2D surface density of the IQMEanh planets at 
r = 0.6, 1.0, 2.5,4. Oay„p (from top to bottom) is displayed. The 
origin of the structure of the standard tp - 1.0 a ./„„ case (second 
from top) was described in iKlev et al] (|2009|) . and we display 
this case here again for better reference. 

The planet located at rp - 0.6ay„p (top panel) is still prone to 
outward migration (see positive torque in Fig.[T]i, but at a slower 
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Fig. 4. Radial torque density r(r) acting on lOMEanh planets 
at rp - 0.6, rp = 1.0, rp - 2.5 and rp - 4.0a jup. The dis- 
tances have been scaled in r-direction to the actual planet loca- 
tion rp.The snapshots of the torque density have been taken in 
the equilibrium state at 80 Jupiter orbits. 



rate. The surface density distribution shows some significant dif- 
ferences compared to the rp = I .Oajup planet. The density in- 
crease ahead and inside of the planet {<p > 180° and r < 0.6) is 
still visible in the rp - 0.6ay„p case, but the density decrease be- 
hind the planet (0 < 180° and r > 0.6) is not that clearly visible. 
Indeed it seems as if the density behind the planet increased in a 
way that the total torque acting on the planet is reduced dramat- 
ically, resulting in reduced positive torque acting on the planet. 

For the planet at rp - 2.5aj„p, where the total torque act- 
ing on the planet is about zero, the density increase ahead of the 
planet is no longer visible in the surface density plot, but the de- 
crease behind the planet is clearly visible. It also seems that the 
planet is able to deplete a larger area around it, possibly owing 
to the onset of gap formation. For even longer distances to the 
central star (e.g. rp - 4.0a jup), the effect becomes ever stronger, 
and the gap is more pronounced. The density increase in front 
of the planet is no longer visible at all. The torque acting on this 
planet is clearly negative, indicating inward migration. With in- 
creasing distance from the host star, the opening angle 6 of the 
spiral arms becomes smaller (see Fig.|5]l, as can be inferred from 
<5 ~ cJx'Kep which scales oc r^O-^^ j-qj. ^^j. temperature gradient 
of T{r) oc r-^'^. 

To analyse the extent of outward migration from our disc 
properties in more detail, it seems useful to compare various 
time scales: the libration, the radiative and the viscous time 
scale in the disc. The necessary unsaturated torques needed 
for sustained outward migration require a pproximately equal 
libration and radiative difiiision times (see iBaruteau & Massed 
(12008 al) : IPaardekooper & Papaloizoul ( l2008h ). For the latter we 
use in our case the rad iative diffusion time scale. Trad- We define 
(iBitsch & Klevll2010l) 
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Fig. 5. Surface density maps for a lOMEanh planet on fixed cir- 
cular orbit in fully radiative discs at four different locations. 
The distance from star to planet changes (from top to bottom): 
rp - 0.6, rp = 1.0, rp = 2.5 and rp = 4.0 (in Jupiter radii). 
Please note the slightly different colour scale for each plot. 
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Fig. 6. Radiative and viscous diffusion time scales that depend 
on the distance from the central star for our standard disc model. 
To compute the time scales we used the density and temperature 
of the midplane at the beginning of the fully radiative simula- 
tions (when the disc is in the r - equilibrium between viscous 
heating and radiative transport/cooling). Libration time scales 
are stated for 20 and 3QMEarih planets. 



For the typical diffusion length s we substitute the vertical disc 
height, i.e. s - H, where we us e H - cJQ. with the sound speed 
Cs- The libration time given by (iBaruteau & Massell2008ah 



T,ih = 87Trp/(3Q.KXs) 



(5) 



where Xs denotes the half-width of the horseshoe-orbit, Xs = 

l.l6rp JqliyjjHIr). Similarly to the radiative diffusion time, 

the viscous time scale is given by Tyi^c = s^ Iv, again with s - H, 
and a constant v. To compute the time scales all required quanti- 
ties are evaluated in the midplane of the unperturbed disc at the 
beginning of the simulations. This applies to the density, temper- 
ature, opacity K{p, T), and the sound speed and Q. 

The three time scales are displayed in Fig. |6l For accretion 
discs that are solely heated internally by viscous dissipation we 
expect in equilibrium Trad ~ Ti/sc- Apparently, this relation is ful- 
filled well (for r < l.5aj„p). We have plotted the libration time 
for two different planet masses, 20 and 3QMEarih- Time scale ar- 
guments suggest a most efficient outward migration for equal 
Tin, and Trad, which is indeed roughly what we find in our 3D 
simulations. However, the overall range of outward migration is 
surprisingly broad. Specifically we find that IQMEanh planets are 
prone to outward migration up to about r x 2.4, where the two 
time scales differ by a factor of 3-4. For 30M£fl,.,/, planets the 
range of outward migration is substantially smaller and centres 
directly around equal libration and radiative time scale. 

We have checked the above estimates of the torques for a sta- 
tionary planet with additional simulations of IQMEanh planets in 
discs, starting at r = l.Qcijup and r - 3.0ajup respectively, which 
were able to move freely inside the disc. The planets gather in 
this case indeed at the zero torque radius (results are not dis- 
played here). 

4. Influence of the disc's mass 

In this section we examine the influence of the disc's mass on 
the migration of low-mass planets embedded in these discs. First 
we compare the relevant physical properties of the discs with 
different masses, and then we investigate the planetary migration 
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in those discs. We then finally discuss convection inside fully 
radiative discs. 



4. 1 . Properties of discs witli different masses 

In our previous work, the disc's mass was fixed to Q.QIMq. We 
now extend the range of disc masses from O.OOSMq to O.O4M0 
(with respect to the standard radial distance, from 0.4 - 2.5). All 
models started locally isothermal with H/r - 0.05 and during 
initial evolution on time this will change to the appropriate equi- 
librium configurations (between viscous heating and radiative 
transport/cooling). 

In Fig. |7] we display the density, temperature, and the as- 
pect ratio of the equilibrium discs for different disc masses at 
r = 1.0. Density and temperature, and H/r are evaluated in the 
disc's midplane. The results are as expected from our previous 
simulations. A higher mass of the disc results in a higher density, 
temperature, and aspect ratio of the disc in the equilibrium state. 
In the isothermal case, a higher aspect ratio of the disc would 
result in a slower i nward migration of a low-m ass planet, see 
iTanaka et alJ (l2002h for linear analysis and e.g. iBitsch & Kiev! 
(1201 Ol) for non-linear simulations. However, low-mass planets in 
fully radiative discs migrate outwards, so that the linear isother- 
mal approach is not valid any more. 

In Fig. [8] the radial distributions of surface density (top) and 
temperature (bottom) are displayed. The temperature has been 
measured in the disc's midplane. By construction, the surface 
density increases for higher disc masses, while it falls off" with in- 
creasing distance to the star, on average according to S(r) oc r"''- 
as expected for a constant v. The surface density profiles for 
the higher mass discs with Mqisc ^ O.O15M0 show some fluc- 
tuations. With increasing disc mass, these fluctuations become 
stronger and reach out to a longer distance from the star While 
they are quite short for Mdisc ^ O.O2M0 and reach only to 
r » 1.3, they become very strong and reach out to r a; 2.3 for 
Moisc = O.O4M0. These fluctuations of the surface density vary 
in time and are related through convective motions in the disc, 
see below. 

The described fluctuations in the surface density can also 
be seen in the temperature profiles of discs with different disc 
masses (bottom panel in Fig.[8]l. The variabilities of the temper- 
ature are not as strong as those for the surface density, neverthe- 
less, they are clearly notable for Mdisc ^ O.OI5M0 and increase 
with the disc's mass. They also, change in time, as does the sur- 
face density. A higher disc mass seems to support stronger fluc- 
tuations that reach farther out into the disc. 

Because a higher disc mass results in a higher aspect ratio 
of the disc, the cut-off of the computed domain (at 7° above 
and below the midplane) might change the structure of the disc. 
Additional simulations with larger did not change the density 
and temperature patterns at 6 - 83° and 9 = 97° for all disc 
masses (see also iKlev et al.lf2009l) . However, a too low bound- 
ary substantially changes the distributions and might therefore 
influence convection in the disc as well. 

The changes in the surface density profiles have a direct 
influence on the migration of an embedded planet. If a planet 
is embedded in a region in the disc where the fluctuations of 
the surface density are very strong (i.e. strong convection), the 
direction of migration might not be clearly determinable, be- 
cause changes in the surface density profile directly influence 
the migration. Stationary gradients in surface density profiles 
can even be used as pla netary traps to collect planetary embryos 
dMorbidelU et al.ll2008h . 
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Fig. 7. Density (top), temperature (middle) and aspect ratio (bot- 
tom) of discs with different masses in the initial equilibrium 
state. All quantities are measured in the disc's midplane at the 
reference distance, r - 1.0. 



4.2. Influences on the migration of low-mass planets 

As mentioned above, a different aspect ratio of the disc will 
change the rate of planetary migration. In the isothermal case, 
a highe r aspect ratio will r esult in slower inward migration in 
theory (Tanaka et al."2002'), which we supported in our previ- 
ous simulations (Bitsch & Kley 2010; Bitsch & Kley 2011). A 
higher disc mass results in a higher aspect ratio and temperature 
of the disc, and we may expect changes in the migration rates. 
When the planet is farther away from the central star, the tem- 
perature and density of the disc are reduced. If the reduction in 
density and temperature is sufficient, the planets stop their out- 
ward migration (see Fig. [TJ. One might now expect that for discs 
with very low masses the torque acting on a IQMEanh planet at 
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Fig. 8. Surface density (top) and temperature (bottom) of discs 
with different masses in the equiUbrium state of the disc. The 
temperature is measured in the disc's midplane. 



r - I .Oajup might become negative. Very high temperatures and 
high densities inside the disc, on the other hand, might influ- 
ence the outward migration as well. Following the formula of 
jtet2010MNRAS.401.1950Pin Equation Q, one might suspect 
stronger torques acting on planets in more massive discs for con- 
stant a, /3 and ^ (increase in surface density overcompensates the 
increase in aspect ratio). 

In the top panel of Fig. |9] the total torque acting on 
IQMEarth plancts on circular orbits embedded in fully radia- 
tive discs with diff'e r ent m asses and the theoretical results from 
iPaardekooper et all (1201 Ol) and iPaardekooper et aP (1201 ih are 
displayed (blue and purple). The torque acting on the planet re- 
mains nearly constant within a small interval around our stan- 
dard disc mass of O.OIMq. For lower disc masses, the torque 
drops off' very rapidly to even negative values for Moisc - 
O.OOSMq, as we expected. For higher disc masses, the torque 
drops off as well. First at a faster rate (to » O.O2OM0), then at a 
slightly slower rate, until it reaches an about zero torque state for 
Moisc - 0.040Mq. This contradicts to our first expectation that 
planets in more massive discs should experience a higher torque, 
the reason may be a change in the temperature gradient and the 
influence of convection. 

When looking at the surface density profile displayed in 
Fig. [8] it is clear that the changes in the surface density may 
disrupt the very sensitive density pattern near the planet. As the 
convection cells in the disc change with time, the torque acting 
on the planet will change as well, giving rise to high fluctua- 
tions/oscillations in the total torque acting on the planet. Hence, 
the torques acting on the planet have been averaged over 20 plan- 



etary orbits. After averaging, the torques acting on planets in 
convective discs show only very low fluctuations. 

F or disc masses around ^ O.OIMq the theoretical formula 
from IPaardekooper et alj ( 1201 Ol) (see eq. [TJ fits our 3D simula- 
tions up to a factor of 25%. For very low disc masses {Mdisc - 
O.OO5M0), however, the fit is not as good. This may be because 
of the reduced disc mass and the consequently changed sur- 
face density distribution (which changes the torque acting on 
the planet), as explained in Section [3] For higher disc masses 
{Mdisc ~ O.OI5M0), the two torque values differ more and more. 
As the torques of our simulations tend to go to zero, the the- 
oretica lly predicted adiabatic torques from IPaardekooper et al.l 
(I2OIOI) become even mo re extended. 

The formula from IPaardekooper et al] (1201 ll) . which in- 
cludes the important effects of viscosity and heat diffusion, dif- 
fers by a factor of three near disc masses of Mdisc ~ 0.0 IM© as 
one could have expected from the res ults presented in Fig . [21 For 
higher disc masses, the torques from IPaardekooper et al.l (1201 ll) 
remain nearly unchanged. The torques from our simulations are 
reduced, but they do not reach negative values. Besides the dif- 
ferences described in Section[3] more differences arise from con- 
vection in the disc, because convection results in fluctuations in 
the torque, which is not considered in the analytical formulae. 
However, convection seems to play a role only for discs with 
Mdisc > O.O2M0. 



4.3. Torque analysis 

In the bottom panel of Fig. |9]the radial torque density r(r) is 
displayed for different disc masses. For the lowest disc mass 
in our simulations, Moisc - O.OOSM©, the usual spike in the 
torque density cannot be seen. The spike in the torque density 
distribution is an indication for a positive torque in a fully ra- 
diative disc (Klevetal. 2009). For higher disc masses (up to 
Moisc ~ O.O2M0), it is is clearly visible. For those disc masses, 
the total torque is indeed positive (see top panel of Fig.|9]l, indi- 
cating outward migration of the embedded protoplanet. 

The torque density for the Moisc = O.O4M0 disc seems to in- 
dicate a total positive torque acting on the planet, and it is indeed 
positive at the moment of the snapshot, but as the fluctuations 
in the surface density change in time, so does the torque act- 
ing on an embedded planet. Therefore the torque density for the 
Moisc = O.O4M0 disc at one single moment during the evolution 
does not necessarily reflect the longterm outcome, if the fluctu- 
ations are to strong. This time variation of the total torque and 
torque density acting on the planet is displayed in Fig. [TO] The 
top panel in Fig. [TOjshows the time evolution of the total torque 
acting on embedded planets for discs with different disc masses. 
For Mdisc < O.OIMq the torque acting on the planet is constant in 
time (after about 50 orbits), while it shows very high fluctuations 
for Mdisc - O.O4M0. In the bottom panel of Fig. [TO] the torque 
density r(r) for a lOMEanh planet in a Mdisc - O.O4M0 disc is 
shown. The torque density is plott ed at different times. Because 
the total torque was fluctuating very much, it is no surprise to 
find these fluctuations for the torque density as well. These fluc- 
tuations, induced by convection, clearly show that a higher disc 
mass disturbs the evolution of the torque. 

In Fig.[TT]the surface density profiles of fully radiative discs 
with disc masses of 0.005, 0.015 and 0.04Mq with embedded 
lOMEanh planets are displayed. The planet in the 0.005Mo disc 
generates a very similar surface density structure compared to 
our standard O.OIM0 disc (second panel in Fig.[5]l. The overall 
density is, of course, reduced (because the disc mass is much 
lower), but the general pattern of the density increases in front 
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Fig. 9. Torque acting on a planet located at 5 AU for differ- 
ent disc masses. Top: Specific total torque acting on planets 
(IQMEarih) embedded in discs with different disc masses. The 
planets embedded in the higher mass discs Molic > O.O2M0 are 
in the convective zone in the disc, so that the torque acting on 
the planet is very noisy and has been averaged over 20 plane- 
tary orbits. Additionally, we over-plotted re sults (blue and pur- 
ple) from the theoretical torque formulae of iPaardekooper et alJ 
(I2010ll201lh for IQMEanh planets. Bottom: Radial torque den- 
sity T{r) acting on the planet for different disc masses. For com- 
parison, the vertical line indicates the location of the maximum 
as found for our standard case. 



of the planet (f> > 180° and r < 1.0 and the decrease behind the 
planet (p < 180° and r > 1 .0 remains intact. However, the density 
structure relevant for outward migration is not as pronounced as 
it should be to result in a positive torque acting on the planet. 

For planets embedded in discs with higher masses, the pic- 
ture is quite different. For a disc mass of O.OISMq (middle pic- 
ture in Fig. fTTT i the density structure in the direction of the star 
(r < 0.9fly„p) is very distorted, but one can still see the density 
increase ahead and the density decrease behind the planet. The 
distortion seems to reduce the torque acting on the planet, but 
the overall torque is still positive, indicating outward migration. 
For an even higher disc mass (bottom picture in Fig. [TT| with 
Moisc = O.O4OM0) the distortions in the disc increase more. The 
density structure, normally seen for low-mass planets in fully 
radiative discs, is no longer visible at all. The distortions are so 
strong that the torque acting on the planet becomes about zero, 
indicating only a low migration rate. 

However, in this last case the torque is in a state of constant 
fluctuations, which complicates realistic predictions about the 
direction of migration in these massive discs. The fluctuations 
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Fig. 10. Top: Specific total torque evolution acting on planets 
(lOMEarih) embedded in discs with different disc masses. The 
planets embedded in the higher mass discs Moisc > O.OIMq are 
in the convective zone in the disc, so that the torque acting on the 
planet is very noisy compared to the low-mass discs. Bottom: 
Radial torque density r(r) acting on the planet embedded in a 
disc with Mdisc = O.O4M0 at different stages of the evolution. 



of the torque have their cause in fluctuations of the density pat- 
terns, which indicates that the convective zone inside the disc is 
enlarged compared to low-massive discs. We observ ed the phe- 
nome non of convection briefly in our previous work dKlev et al.l 
120091) for discs with Moisc - O.OIMq as well, but the convective 
zone did not reach the planet, and thus did not disturb the density 
pattern around the planet. 

In Fig. [12] we display the radiative diffusion time scale and 
the libration time scale for discs with different masses. In or- 
der to keep the torques acting on the planet unsaturated (to 
evoke outward migration), the libration an d radiative diffusion 
time scale ne ed to be approximately equ al jBaruteau & Massea 
2008 a; Paarde kooper & Papaloizoul2008h . However, for convec- 
tive discs the equilibrium state through cooling through convec- 
tion is reached only when Trad > t"i7.k as observed in Fig. [12] The 
time scales for the 0.005 Mq disc indicate outward migration in 
a region at r » l.Oajup, but at r = l.Oay,,,, the time scales differ 
by a factor of 4, indicating that the torques are not kept unsat- 
urated in this region of the disc, which in turn indicates inward 
migration (as presented in in the top figure in Fig.|9]l- 

For the O.OISMq disc the time scales are nearly identical at 
r - l.Oajiip, which indicates outward migration (as can be seen 
in the top of Fig. |9]l- However, for longer distances to the central 
star, the time scales start to differ, which indicates inward mi- 
gration. In the O.O4OM0 disc, the time scales differ by a factor 
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Fig. 11. Surface density maps for planets on fixed circular or- 
bits in fully radiative discs with different disc masses (from top 
to bottom); Mdisc - O.OOSMq, Mdisc = O.OISMq and Mdisc - 
0.040Mo. The disruption in the surface density patterns for the 
higher mass discs are caused by convection inside the disc. 



of 3 at r = 1 -Oajup, which indicates inward migration. However, 
the measured torque acting on the planet is positive, indicating 
slow outward migration. But because the planet is embedded in 
a highly convective region in the disc, it is very difficult to pre- 
dict the motion of the planet correctly by considering only the 
time scales. 
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Fig. 12. Time scales depending on the distance from the cen- 
tral star for O.OOSMq, O.OIOMq, O.OlSMo and O.O4OM0 discs. 
To compute the time scales we used the density and temperature 
of the midplane at the beginning of the fully radiative simula- 
tions (when the disc is in the r - 6 equilibrium between viscous 
heating and radiative transport/cooling). 



4.4. Orbital evolution 

In Fig. [13] the evolution of semi-major axis for IQMEanh plan- 
ets in isothermal and fully radiative discs with different disc 
masses is displayed. The isothermal reference simulations are 
started from a Hjr = 0.037 disc, which represents the Hjr value 
at the planets starting location in the fully radiative regime. In 
the isothermal disc, no convection is present and therefore the 
embedded planet migrates as expected. A higher disc mass re- 
sults in a faster inward migration. However, it seems that for the 
Mdisc - O.O4M0 disc the type-III-migration regime is hit, be- 
cause the planet moves inwards very fast. 

For a planet on a fixed circular orbit in a fully radiative 
disc with Mdisc - O.O4M0 we determined a positive torque (see 
Fig. |9]l by averaging in time. However, the total torque acting 
on the planet undergoes strong fluctuations in time (see Fig.fTOli. 
When embedding a IQMEarth planet in such a highly convec- 
tive disc, the evolution pattern should to some extend reflect the 
strong fluctuations in the torque acting on a planet on a fixed 
orbit. Indeed, the planet experiences some small kicks in its evo- 
lution pattern (see Fig. [T3T l. Interestingly, the planet migrates 
inwards despite the positive torque acting on a fixed planet. It 
seems that in the convective region of high-mass discs, the mea- 
surement of the torque for planets on fixed orbits is not as re- 
liable as for planets in low-mass discs. Additionally, the time 
averaging may have to be performed over a longer time span. 

4.5. Convective zone 

In order to investigate whether convection is actually present in 
the disc, we display the velocities in z-direction (out of mid- 
plane of the disc) for difl'erent disc masses {Moisc - 0.005 Mq, 
Moisc - 0.015Mq and Moisc - O.O4M0) in the disc's midplane 
in Fig. [14] These plots represent simulations with a whole disc, 
meaning 83° < 6 < 91°. As the surface density plots indi- 
cated, there are no disrupted areas in the velocity patters for the 
Moisc = O.OO5M0 simulations. Therefore for this low disc mass 
we observe no convection near the planet. 

However, for the Moisc = O.OI5M0 case the surface density 
patterns already indicated that convection is possible in the disc 
inside of the planet's distance to the central star The velocity dis- 
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Fig. 13. Evolution of semi-major axis for 20MEarth planets (start- 
ing at r = l.Oajiip) in isothermal and fully radiative discs with 
different disc masses. 



tributions confirm this result, lines with positive velocities (indi- 
cating a flow towards the upper boundary of the disc) alternate 
with lines with negative velocities (indicating a flow towards the 
midplane of the disc). These flows are a clear indicator of mo- 
tion caused by convection in the disc. In the Moisc = O.O4OM0 
case, the fluctuations in the surface density increased dramat- 
ically, as did the changes in the velocity pattern. Alternations 
between positive and negative velocities have increased and in- 
dicate a very strong convective region that disturbs the torque 
acting on the planet (and therefore it's migration). The flow pat- 
tern in the disc is very erratic, making it absolutely necessary 
to average the torque acting on an embedded planet for many 
orbits. 

For simulations that cover only the upper half of the disc, 
the convection cells inside the disc end at the disc's midplane, 
but in reality these convective cells continue to the lower half of 
the disc. However, that convection inside the disc changes the 
behaviour of the planet disc interactions can also be seen for 
simulations containing only one half of the disc. 

In Fig. [15] the velocity in z direction is displayed for a half- 
size disc (only the upper half of the disc is computed, top figure) 
and for a full disc. These computations have been performed 
only for fully radiative axisymmetric 2D discs (in r-6 direc- 
tion) with a disc mass of Mdisc - O.O4M0 without an embedded 
planet. The convection in the disc is clearly visible. In both sim- 
ulations the convection cells in the disc become more symmetric 
for distances longer than r > 1 .25fly„p, meaning that the veloc- 
ity changes from positive to negative only in radial direction and 
not in the vertical direction as well. For shorter distances to the 
central star the convection cells are very irregular in both cases. 

When putting a lOMEanh planet in the 0.04Mo discs, the 
structure of the convective cells changes as the disc gets dis- 
turbed by the planet. Without the planet, the convective cells 
are very regular for distances to the central star exceeding 
r > l.25ajup. With an embedded planet these regular struc- 
tures break down and become very irregular, as can be seen in 
Fig. [161 This effect may caused by to the wake created by the 
planet inside the disc, which acts as an additional heat source. 
At shorter distances to the central star, the structure is irregular 
with or without an embedded planet. 

The velocity patterns for the one sided and two sided discs 
show a little difference. It seems that the fluctuations in veloc- 
ity are more centred in midplane for the one sided disc, while 
they seem to be located near the upper and lower boundaries for 
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Fig. 14. Vertical velocities in z-direction for planets on fixed cir- 
cular orbits in fully radiative discs with different disc masses 
(from top to bottom) in the disc's midplane; Moisc - O.QOSMq, 
Moisc - O.O15M0 and Mohc - O.O4M0. The disruption in the 
velocity patterns for the higher mass discs are caused by convec- 
tion inside the disc. A positive velocity simulates outward flows, 
while a negative velocity indicates a flow towards the disc's mid- 
plane. 



the two sided discs. This effect has several reasons. In the one 
sided disc, material is reflected at the midplane of the disc, which 
might lead to an increase of fluctuations near the midplane. In the 
two sided disc, material can flow through the midplane, so that 
the fluctuations near the midplane are reduced. Furthermore, the 
one sided disc might be unrealistic if convection is present in the 
disc. 

However, the general structure changes when both halves 
of the disc are computed, independent of an embedded planet. 
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Fig. 15. Velocities in z-direction for 0.04Mo discs without 
planet. In the top panel only the upper half of the disc was sim- 
ulated, while in the bottom panel both sides of the disc were 
simulated (with twice the number of grid cells in direction). 
The other simulation parameters are identical. 




Fig. 16. Velocities in z-direction for O.O4M0 discs with embed- 
ded lOMEarth planet. In the top panel only the upper half of the 
disc was simulated, while in the bottom panel both sides of the 
disc were simulated (with twice the number of grid cells in 9 
direction). The other simulation parameters are identical. 



The convection cells are now moving through the midplane of 
the disc, which was not possible for sim u lations of only on e 
half of the disc, see also lKlev et aP (Il993h : iKIahr et al.l (Il999h . 
Therefore, the surface density structure in the two sided disc 
(with Miiisc = O.O4M0) is slightly different compared to the one 
sided disc. In the two sided case, the fluctuations in the surface 
density continue only to a; 1 .5fly„p, while they covered the whole 
disc in the one sided case. The temperature profiles, on the other 
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Fig. 17. Specific torque acting on 2QMEarrh planets embedded in 
0.02, 0.025, 0.03 and O.O4M0 discs at the following distances 
rp = 2.0, 2.5, 3.0 and 4.0a y„p, respectively. 



hand, show no difference at all. We therefore only expect little 
change in the torque acting on planets embedded in one or two 
sided discs at rp - I .Oajup. 

Simulations of embedded lOMEanh planets in fully radiative 
discs that cover both sides of the disc show only very small dif- 
ferences in the velocity pattern in mid-plane of the disc. The 
fluctuations in time of the torque acting on the planet embed- 
ded in Mdisc - 0.04Mq and M^usc = 0.03Mo discs are compa- 
rable (with only small differences in the amplitude of the fluc- 
tuations) for both simulations, confirming our previous assump- 
tions. Simulations of planets embedded in lower mass discs show 
no difference at all (simulations not displayed here), because the 
convective region does not reach to the planet at all. 

Considering that a longer distance from the central star re- 
sulted in a turn from outward to inward migration for a lOMEarth 
planet (see Fig. [TJ because of a reduction in temperature and 
density at the given location of the planet, one might argue that 
the zero-torque distance from the central star might be at longer 
distances for higher disc masses. Moreover, as the disrupting 
convective zone in massive discs reaches farther out from the 
star, outward migration might be possible at larger radii in more 
massive discs, because the disc's convective zone stops at longer 
distance to the central star and the density is still high enough to 
produce the surface density patter needed for outward migration. 

Additional simulations with lOMEarth planets in 0.02, 0.025, 



0.03 and 0.04Mo discs with r 



p 



2.0, 2.5, 3.0 and 4.0fl 



Jap: 



re- 



spectively, confirm our assumption (see Fig. [TTl i. The torques 
acting on those planets are positive, indicating outward migra- 
tion, and show no fluctuations in time. The surface density plots 
also show no sign of convection in the disc at the location of 
the planet (not displayed here). It seems that outward migration 
is therefore possible to farther distances from the central star in 
more massive discs. 

The picture of convection in our disc would change when in- 
cluding stellar irradiation because it would heat the surfaces of 
the disc in contrast to the applied cooling right now. This would 
result in less convection in the disc. Because the convective re- 
gion is a result of the higher surface density (increasing Trad) 
and viscosity in the disc, a reduction in viscosity could prevent 
convection in the disc. However, a reduction of viscosity also 
reduces the torque of an embedded planet, so that outward mi- 
gration might not be possible any more, even for low-mass discs. 
The influence of viscosity will also be addressed in a next paper 
in much more detail. 
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In self-gravitating discs, the torque acting on an embed- 
ded planet can differ by a facto r of two co mpared to non 
self-gravitating discs, as shown bv lBaruteau & M asset (2008b). 
These authors also state that self-gravity has no effect on the 
corotation torque in the linear regime, but our 3D simulations are 
in the non-linear regime. Therefore the influence of self-gravity 
on planet migration in fully radiative discs should be investigated 
in the future. 

The Toomre stability criterion c an be used to estimate the 
stability of discs against self-gravity (iToomrdl 19641) . The stabil- 
ity parameter reads 



G = 



7TI.G 



(6) 



where Cs is the sound speed in the disc, K^p is the epicyclic fre- 
quency, which for Keplerian discs is approximately equal to the 
angular frequency Q, E is the surface mass density and G is the 
gravitational constant. In order to achieve stability in discs, the 
stability parameter must he Q ^> 1 . For all disc masses used in 
this work, this criterion is fulfilled well, so that the discs are not 
gravitationally unstable. 

Because convection is a 3D effect, 2D simulations (in r-0 
direction in the midplane) of fully radiative discs with high discs 
masses (M,/,,^ > O.O2M0) cannot capture this effect. Therefore 
planets embedded in these 2D simulations will not be exposed 
to these fluctuations and might therefore be inaccurate near the 
central star because of convection in the disc. 



5. Summary and conclusions 

We performed full 3D radiation hydrodynamical simulations of 
low-mass planets embedded in accretion discs at different dis- 
tances to the central star and for various disc masses. 

In the first sequence of our simulations we changed the plan- 
etary distance to the central star of embedded planets on circular 
orbits. With increasing distance to the central star, the torque 
acting on IQMEanh planets embedded in fully radiative discs be- 
comes even more reduced and it reaches negative torques for 
longer distances. We find an equilibrium, zero-torque distance, 
to the central star for IQMEarth planets at r a; 2.4a y,,,,. This equi- 
librium distance varies with the planetary mass (for ISMEanh 
planets it is r ^ 1 .9aj„p and r ^ 1 Aajup for 'iQMEarth planets), 
indicating that a quite extended region in the disc might act as a 
feeding zone to create even larger planetary cores. The concept 
of equilibrium radius (zero torque radius) f or planetary embry os 
in fully unsaturated discs has been stated in lLvra et al.l ( 1201 Ob as 
well and it could easily act as a feeding or collection zone for 
planetary embryos. 

Planets embedded in fully radiative discs migrating outwards 
create a very sensitive pattern in the surface density distribution. 
Ahead and inside of the planet a density increase is visible (see 
second panel from the top in Fig. |5]l, which shrinks with in- 
creasing distance to the central star. This density enhancement 
is indeed accountable for the positive torque acting on the planet 
(also visible as a spike in the radial torque density distribution in 
Fig. |4]i, but as the distance to the central star increases, this ef- 
fect becomes less, so that it cannot overcompensate the negative 
Lindblad torques any more, which results in inward migration. 

We compared our results to the recently developed 
I 1 I 1 I 1 

torque formulae by Paardekooper et al. (2010, 2011) and 
iMasset & Casolil (2010). P aardekooperet a l. (2010) includes 
just the fully unsaturated torques in the inviscid and adiabatic 
case, shows no torque reversal option for constant /? and is as 



such unphysical when comparing it with the long-term evolu- 
tion of planets. However, as in our simulations /? changes in the 
disc with distance to the central star (fiir - 1.0fly„p) = 1.7 to 
f5{r - S.Oajiip) - 0.33), the torque given by the formula be- 
comes negative for large distances to the central star. But, overall 
the match of the torque of the formula with our simulations is not 
good. This formula is only valid in the first orbits after the planet 
is embedded when the torques are still unsaturated, however, the 
torques do saturate in time. 

However, as expect ed, the improved version of 
iPaardekooper et all (1201 ll) that includes viscous and heat 
diffusion describes our results more accurately. It also features 
a transition from positive to negative torques, but the negative 
torques at large distances to the central star are about a factor 
of two larger than in our simulations. Even though an exact 
match of the torques has not been achieved, the trend seems to 
be caputured quite well by the formula. 

The torque formulae are derived for 2D discs with thermal 
diffusion operating only in the disk's midplane, while our sim- 
ulations are 3D, and take full account of vertical diffusion as 
well, which can change the structure of the disc. The formulae 
were also derived for given gradients in temperature and sur- 
face density, but in a real disc the temperature and surface den- 
sity profiles are disturbed when a planet is embedded in a disc. 
Because the formulae were derived and checked for a SMEanh 
planet (with about 20 per cent agreement), the disturbances of 
such a small planet in a disc are much weaker than for our em- 
bedded IQMEanh planet, which may give rise to more significant, 
non-linear disturbances in the temperature and density profiles. 
All of this may lead to differences between our simulations and 
the theoretical for mulae in the torque acting on the planet. 

The formula of iMasset & Casolil (1201 Ol) does not match our 
simulations that well. The torques from the formula are al- 
ways negative, which is in contrast to our simulations, where 
we find outward m igration for r < 2 . 5a./»p . For long distances 
to the central star, IMasset & Casolil (1201 Oh match better than 
IPaardekooper et aP (1201 1|). as the Lindbl a d torq ue is reduced. 
The Lindblad torque of IMasset & Casolil (12010 ') also matches 
quite well with the isothermal torque of iD'A ngelo & Lubo^ 
(I2OIOI) . when taking the factor y into account due to the dif- 
ferent sound speeds in isothermal and fully radiative discs. The 
overall trend to h ave torques from positive to negative seems 
best achieved with'Paa rdekooper et al.l (201 1"). but a quantitative 
agreement is still lacking. In Appendix lAl we discuss the influ- 
ence of the smoothing length on the formulae. 

For increasing disc masses, the temperature, density (in mid- 
plane), and aspect ratio of the disc increases in the equilibrium 
state where viscous heating and radiative transport/cooling are 
in balance. The convective zone in the inner discs stretches far- 
ther out with increasing disc mass, resulting in high fluctuations 
of the surface density in our computed domain for discs with a 
mass higher than Mdisc ~ 0.02Mo. 

Starting from a Mdisc - O.OIM© disc, the torque acting 
on embedded IQMEanh planets decreases for increasing disc 
masses. As the disc mass increases, the convective zone in the 
disc stretches farther out from the central star and influences 
planetary migration. The fluctuations in the disc's density disrupt 
the torque acting on the planet on a stationary orbit for high-mass 
discs in a way that the torque is very irregular and shows high 
fluctuations as well, making it difficult to determine the correct 
direction of migration. For lower disc masses, the torque reduces 
as well, presumably because of the same reasons as the torque 
is reduced for longer distances to the central star in a discs with 
Md,„ = O.OlMo. 
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The formula in iPaardekooper et al.l (1201 ih fits within a fac- 
tor of three with our 3D simulations for planets in discs with 
Mciics ~ O.OIM0. For higher disc masses, the difference between 
the formula and our 3D simulations becomes smaller However, 
the 3D simulatio ns show a decreasing torque for increasing disc 
mass, while the IPaardekooper et alT (1201 ll) formula increases 
slightly. As the disc mass increases, viscous heating increases 
and cooling becomes inefficient, which results in a structure sim- 
ilar to an adiabatic disc. In adiabatic discs the corotation torques 
saturates, resulting in a lower torque acting on the planet, hence 
the drop of the torque for increasing disc masses. Convection 
certai nly plays a role in more ma ssive discs, but it is unaccounted 
for in lPaardekooper et al.l (1201 ll) . In more massive discs the con- 
vective zone reaches longer distances from the central star, dis- 
rupting the density pattern near the embedded planet and thus 
creating fluctuations in the torque acting on the planet. These 
disruptions in the density pattern are caused by the convective 
cells evolving in the disc. These cells also change in time, giving 
rise to the stronger fluctuations of the torque. 

Convection is ineffic ient for transporting angular momentum 
dKIev et al.ll 19931: iLesur & Ogilvie 2010), but the influences of 
convection on planetary migration are very dramatic, because 
a planet close to the star in the convective zone of the disc is 
essentially disrupted. The direction of migration is not clearly 
determinable any more, but when the planet is farther out in 
the massive disc and the convection fades away, the direction 
of migration is easy to specify, indicating outward migration. 
Therefore the zero-torque radius for migration lies farther out in 
more massive discs. 

Convection is also a 3D effect only and cannot be simulated 
in 2D discs (in r-(p direction in the midplane). Two-dimensional 
simulations of planets in massive discs {Mdisc > O.O2M0) might 
therefore be inaccurate near the central star, because the effects 
of convection are not considered. 



Ap pendix A: Comparison witii IPaardekooper et al 

(l201lh 



In Section (O we compared our numerical results to the analyti- 
cal torque formula derived in lPaardekooper et alj (1201 lb . As the 
derivation is rather cumbersome, we present here a brief sum- 
mary of the relevant contributions to the total torque acting on 
a planet, so that the reade r can follow our ca l culatio ns. In our 
notation we closely follow IPaardekooper et al.l (1201 ll) . 

The total torque acting on a low-mass planet consists of two 
main contributions the Lindblad torque, Ti, plus the corotation 
torque, Fc 



Ffo, — Tr + Fr 



(A.l) 



The Lindblad torque is caused by the action of the induc ed spiral 
arms and is given as (IPaardekooper & Papaloizoull2008h 



yFi/Fo^ -2.5-1.7/3 + 0.10'. 



(A.2) 



where a denotes the negative slope of the surface density profile 
S oc r^", fi refers to the slope of the temperature profile T oc r^P, 
and y is the adiabatic index of the gas. 

It is important to notice that all torques listed here are nor- 
malized to 



(ir^-; 



■;"j 



with q the planet/star mass ratio, Ep the surface density at the 
planet's location, Q.p the rotation frequency of the planet in the 
disc, and h is in this appendix the isothermal relative disk height. 



The corotation torque is now split into the barotropic part 
and an entropy-related part: 



F,. = F. 



,bai'o 



+ Tr 



where the first part applies to barotropic flows where the pressure 
only depends on the density, and it depends on the gradient of the 
vorticity in the flow; the second part relates to the variations of 
entropy. Each of them is split again into a linear contribution and 
a so-called horseshoe drag contribution. This separation is nec- 
essary because the two parts are affected differently by the diffu- 
sion processes. The barotropic part of the (non-Unear) horseshoe 
drag is given by 

yTh,,barolU^l.\{l.5-a) (A.3) 

and the entropy -related part of the horseshoe drag is given by 



y^hs.enll'^O =7.9- 



r 



(A.4) 



where ^ - p - (y - 1 .Q)a is the negative of the power-law index 
of the entropy. We note th at the total torque formula given by 
IPaardekooper et all (l2010t) . as summarized in Eq. ([1]), is exactly 
the sum T,ot -Ti+ Ths.baw + ^hs.em- 

The barotropic part of the linear corotation torque reads as 



y^cMnJmrol^a = 0.7(1.5 - a). 



(A.5) 



and the entropy-related part of the linear corotation torque is 
given by 



y^cjm,entl^0 = (2.2 — )^. 



(A.6) 



Owing to the difference between the isothermal and adiabatic 
sound speed, differences in the torque arise. To compensate for 
this, the adiabatic index y should be replaced by an effective y: 



2Qy 



yQ + \^2 ^Jif-Q^ + 1)2 - i6e-(r^ + 27^22 Z^ 



yeff^ 



so that all y's in the previous equations ( IA.2l to lAT6i l have to be 
replaced by y^//. The parameter Q is given by 



e = 



2xpQ.p 



2Xf 



2>hcl 3h^rlD.p ' 
where h-Hjr and 



Xp 



16y(y-l)trpr^ 
3KpplH-pQo 



with K beeing the opacity and cr the Stefan-Boltzm ann con- 
stant. Please note that in IPaardekooper et all (1201 ll) a factor 
of 4 is missing. The final correction relates to the non-ideal 
effects of viscosity and heat transfer, which both have to be 
present to avoid the saturation of the corotation torque. The 
barotropic part of the horseshoe drag is not affected by thermal 
diffiision and is only deter mined by the viscosity. According to 
IPaardekooper et all (1201 ih it can be written as 

^Chnro = ^hs.baroF(Pv)G(pv) + (1 - K{py))Yc,Un.baro, 

where F;,5j,„„, and Tcjm,haro are given by equations ( IA.3l and [A.5b . 
but now with y — » yeff- F{p) (eq. IA.7l i governs saturation and 
G{p) (eq. IA.8I 1 and K{p) (eq. |A.9t govern the cut-off at high 
viscosity. 
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For the non-barotrop ic, entropy-related corotation torque 
iPaardekooper et alJ (1201 ih find 



A c,e, 



= ^hs.entF{pv)F(p^) y/G(py)G(pj^) 

+ ^J{l-K{pM^-K(p,W,,n„,,„, , 



where r/„,„« and Tc^im,ent are given by equations IA.4I and IA.6I 
again with y — » Jeff, and p^ is the saturation parameter associ- 
ated with thermal diff'usion. 

The function F{p) is given by 



F{p)^ 



1 



l+(p/ 1.3)2- 
The function G{p) is given by 



G{p) 



m^-^f^p^i'- 



25 V 8 



\4/3 






for p < 
for p > 



A5n 
45Jr 



(A.7) 



(A.8) 



The function K{p) is given by 



K{p)^ 



i6/452\3/4 3/2 
25 \ 28 / '^ 

_ 9 f' 28 \"^^-^„-8/3 
15\A5itI i' 



for /^<VS 



for 



(A.9) 



P^ \^4^ 



The parameters pv and p^, relate to the strength of viscosity 
and thermal difFusivity, and are given by 



Pv = 3 



Px 



r-pQpxl 



27rv„ 



'^"Xp 



where Vp is the kinematic viscosity and Xp the thermal conduc- 
tivity at the planet location, x, is the half width of the horseshoe, 
given by 



- 1-1 /0-4\ 

'<!ff 



1/4 



The scaling with e//i breaks down for small softening (e//i < 
0.3). All these contributions have to be substituted into equation 
dA.lb to calculate the total torque acting on the planet. In deriv- 
ing these formulae one has to use a description for the smoothing 
of the gravitational potential. Here, a standard e potential has 
been assumed, and a smoothing length of ejh - 0.4 has been 
made. 

To compare with our simulations, we have used the follow- 
ing parameter in the formulae above a - 0.5, /3 = 1.7,^ = 1.5. 
To evaluate the viscosity parameter py we use v = 10"^. Please 
note t hat v,, ^ Xo for di scs in radiative equilibrium. 

In lKlev et alJ (l2009l) we pointed out that a different planetary 
smoothing results in different torques. This phenomenon can be 
up to a factor of two for the e potential with rs„, = 0.5 and the 
cubic potential with rs„, - 0.5. In Section (O we compared our 
numerical simulations to the smoothing with e/h = 0.4. 

Because our cubic potential with r,,,, is deeper in the vicinity 
of the planet, we used e/h - 0.3 , e/h = 0.35 and e/h - 0.5 
for the IPaardekooper et al.l (1201 ih formula as well. The results 
are shown in Fig. lA.ll In the e/h = 0.5 case, the torque is nega- 
tive for all planetary distances, which is not a match at all for our 
simulations, but interestingly agrees better with Masset & Casolil 
(I2OIOI) . For e/h - 0.3, the torque is much higher in the positive 
torque regime and much lower in the negative torque regime. 



0.00015 r 
0.0001 [ „ M 
5e-05 

-5e-05 - 
-0.0001 - 

-0-00015 - 3D Sim, 20ME,rth 

eps/h=0.4 ■ ■ X- ■ 
-00002 - eps/h=0.3 « 

eps/h=0.5 --B-- 

-0.00025 - eps/h=0.35 • •■ • 

-0.0003 I ' ' ' ' ' ' — 

0.5 1 1.5 2 2.5 3 3.5 

distance to star [ajup] 




4.5 



Fig.A.l. Specific torque acting on lOMEanh planets em- 
bedded in a 0.0 IMp di sc. Overplotted are the results of 
IPaardekooper et aP (1201 ll) using various smoothings lengths 
from e/h = 0.3 and e/h - 0.5. 



For the e/h = 0.35 case, the formula matches our simulations 
quite well in the positive torque regime (r < 2.5a j„p). However, 
for larger distances to the central star, the torque of the formula 
is about a factor of two to four larger than the torque from our 
simulations. Nevertheless, the match for r < 2.5aj,ip is quite 
good. Obviously, the smoothing of the planetary potential seems 
to have a huge effect on the torque acting on the planet. The ef- 
fect seems much stronger in the formula compared to different 
smoothings in our 3D hydrodynamical simulations dKlev et al.l 
I2OO9I) but they are based on 2D simulations, where a larger im- 
pact is to be expected. 
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